From 42d98a1543b0a2ee432610b0554abc134e4c4a57 Mon Sep 17 00:00:00 2001 From: robertl Date: Wed, 17 Jun 2009 01:29:14 +0000 Subject: [PATCH] Amit Gil adds support for Israeli Transverse Mercator/New Israeli Grid for Naviguide. --- jeeps/gpsdatum.h | 8 +- jeeps/gpsmath.c | 379 +++++++++++++++++++++++++++++++--- reference/route/naviguide.gpx | 40 ++-- 3 files changed, 376 insertions(+), 51 deletions(-) diff --git a/jeeps/gpsdatum.h b/jeeps/gpsdatum.h index 7cbe23ec1..60f14619a 100644 --- a/jeeps/gpsdatum.h +++ b/jeeps/gpsdatum.h @@ -44,8 +44,7 @@ GPS_OEllipse GPS_Ellipse[]= { "WGS66", 6378145.000, 298.25 }, { "WGS72", 6378135.000, 298.26 }, { "WGS84", 6378137.000, 298.257223563 }, - { "Clarke 1880(Benoit)", 6378300.789,293.4663155389811 }, - + { "Clarke 1880 (Benoit)", 6378300.789, 293.466 }, }; @@ -185,7 +184,8 @@ GPS_ODatum GPS_Datum[]= /* 121 */ { "Sweden", 4, 424.3, -80.5, 613.1 }, /* 122 */ { "GDA 94", 21, 0, 0, 0 }, /* 123 */ { "CH-1903", 4, 674, 15, 405 }, -/*124 */ { "Palestine 1923", 27, -235, -85, 264}, +/* 124 */ { "Palestine 1923", 27, -235, -85, 264 }, +/* 125 */ { "ITM (Israeli New)", 21, -48, 55, -52 }, { NULL, 0, 0, 0, 0 } }; @@ -224,6 +224,8 @@ GPS_ODatum_Alias GPS_DatumAlias[] = { "WGS-72", 117 }, { "WGS84", 118 }, { "WGS-84", 118 }, + { "Israeli", 124 }, + { "D_Israel_new", 125 }, { NULL, -1 } }; diff --git a/jeeps/gpsmath.c b/jeeps/gpsmath.c index e0c31e118..6a5ee3df2 100644 --- a/jeeps/gpsmath.c +++ b/jeeps/gpsmath.c @@ -733,6 +733,329 @@ void GPS_Math_Swiss_EN_To_WGS84(double E, double N, double *lat, double *lon) } +/* @func GPS_Math_Cassini_LatLon_To_EN ********************************** +** +** Convert latitude and longitude to Cassini transverse cylindrical projection +** easting and northing +** +** @param [r] phi [double] latitude (deg) +** @param [r] lambda [double] longitude (deg) +** @param [w] E [double *] easting (metre) +** @param [w] N [double *] northing (metre) +** @param [r] phi0 [double] latitude of origin (deg) +** @param [r] M0 [double] central meridian (deg) +** @param [r] E0 [double] false easting +** @param [r] N0 [double] false northing +** @param [r] a [double] semi-major axis +** @param [r] b [double] semi-minor axis +** +** @return [void] +************************************************************************/ +void GPS_Math_Cassini_LatLon_To_EN(double phi, double lambda, double *E, + double *N, double phi0, double M0, + double E0, double N0, double a, double b) +{ + double p2; + double po2; + double a2; + double b2; + double e2; + double e4; + double e6; + double AM0; + double c0; + double c1; + double c2; + double c3; + double om0; + double A0; + double A1; + double A2; + double A3; + double j; + double te4; + double phi0s2; + double phi0s4; + double phi0s6; + double lat; + double x; + double E1; + double E2; + double E3; + double E4; + + double phis; + double phic; + double phit; + double phis2; + double phis4; + double phis6; + double RD; + double dlam; + double NN; + double TT; + double WW; + double WW2; + double WW3; + double WW4; + double WW5; + double CC; + double MM; + + + lambda = GPS_Math_Deg_To_Rad(lambda); + phi0 = GPS_Math_Deg_To_Rad(phi0); + phi = GPS_Math_Deg_To_Rad(phi); + M0 = GPS_Math_Deg_To_Rad(M0); + + + p2 = (double)GPS_PI * (double)2.; + po2 = (double)GPS_PI / (double)2.; + + a2 = a*a; + b2 = b*b; + e2 = (a2-b2)/a2; + e4 = e2*e2; + e6 = e2*e4; + + te4 = (double)3.0 * e4; + j = (double)45. * e6 / (double)1024.; + c0 = (double)1.0-e2/(double)4.-te4/(double)64.-(double)5.*e6/(double)256.; + c1 = (double)3.*e2/(double)8.+te4/(double)32.+j; + c2 = (double)15.*e4/(double)256.+j; + c3 = (double)35.*e6/(double)3072.; + + lat = c0*phi0; + phi0s2 = c1 * sin((double)2.*phi0); + phi0s4 = c2 * sin((double)4.*phi0); + phi0s6 = c3 * sin((double)6.*phi0); + AM0 = a * (lat-phi0s2+phi0s4-phi0s6); + + om0 = (double)1.0 - e2; + x = pow(om0,(double)0.5); + E1 = ((double)1.0 - x) / ((double)1.0 + x); + E2 = E1*E1; + E3 = E1*E2; + E4 = E1*E3; + A0 = (double)3.*E1/(double)2.-(double)27.*E3/(double)32.; + A1 = (double)21.*E2/(double)16.-(double)55.*E4/(double)32.; + A2 = (double)151.*E3/(double)96.; + A3 = (double)1097.*E4/(double)512.; + + + dlam = lambda - M0; + if(dlam>GPS_PI) + dlam -= p2; + if(dlam<-GPS_PI) + dlam += p2; + + phis = sin(phi); + phic = cos(phi); + phit = tan(phi); + RD = pow((double)1.-e2*phis*phis,(double).5); + NN = a/RD; + TT = phit*phit; + WW = dlam*phic; + WW2 = WW*WW; + WW3 = WW*WW2; + WW4 = WW*WW3; + WW5 = WW*WW4; + CC = e2*phic*phic/om0; + lat = c0*phi; + phis2 = c1 * sin((double)2.*phi); + phis4 = c2 * sin((double)4.*phi); + phis6 = c3 * sin((double)6.*phi); + MM = a * (lat-phis2+phis4-phis6); + + *E = NN*(WW-(TT*WW3/(double)6.)-((double)8.-TT+(double)8.*CC)* + (TT*WW5/(double)120.)) + E0; + *N = MM-AM0+NN*phit*((WW2/(double)2.)+((double)5.-TT+(double)6.*CC) * + WW4/(double)24.) + N0; + return; +} + + + + +/* @func GPS_Math_Cassini_EN_To_LatLon ********************************** +** +** Convert Cassini transverse cylindrical easting and northing projection +** to latitude and longitude +** +** @param [r] E [double] easting (metre) +** @param [r] N [double] northing (metre) +** @param [w] phi [double *] latitude (deg) +** @param [w] lambda [double *] longitude (deg) +** @param [r] phi0 [double] latitude of origin (deg) +** @param [r] M0 [double] central meridian (deg) +** @param [r] E0 [double] false easting +** @param [r] N0 [double] false northing +** @param [r] a [double] semi-major axis +** @param [r] b [double] semi-minor axis +** +** @return [void] +************************************************************************/ +void GPS_Math_Cassini_EN_To_LatLon(double E, double N, double *phi, + double *lambda, double phi0, double M0, + double E0, double N0, double a, double b) + +{ + double p2; + double po2; + double a2; + double b2; + double e2; + double e4; + double e6; + double AM0; + double c0; + double c1; + double c2; + double c3; + double om0; + double A0; + double A1; + double A2; + double A3; + double j; + double te4; + double phi0s2; + double phi0s4; + double phi0s6; + double lat; + double x; + double E1; + double E2; + double E3; + double E4; + + double dx; + double dy; + double mu; + double mus2; + double mus4; + double mus6; + double mus8; + double M1; + double phi1; + double phi1s; + double phi1c; + double phi1t; + double T; + double T1; + double N1; + double R1; + double RD; + double DD; + double D2; + double D3; + double D4; + double D5; + double tol; + + M0 = GPS_Math_Deg_To_Rad(M0); + phi0 = GPS_Math_Deg_To_Rad(phi0); + + p2 = (double)GPS_PI * (double)2.; + po2 = (double)GPS_PI / (double)2.; + + a2 = a*a; + b2 = b*b; + e2 = (a2-b2)/a2; + e4 = e2*e2; + e6 = e2*e4; + + te4 = (double)3.0 * e4; + j = (double)45. * e6 / (double)1024.; + c0 = (double)1.0-e2/(double)4.-te4/(double)64.-(double)5.*e6/(double)256.; + c1 = (double)3.*e2/(double)8.+te4/(double)32.+j; + c2 = (double)15.*e4/(double)256.+j; + c3 = (double)35.*e6/(double)3072.; + + lat = c0*phi0; + phi0s2 = c1 * sin((double)2.*phi0); + phi0s4 = c2 * sin((double)4.*phi0); + phi0s6 = c3 * sin((double)6.*phi0); + AM0 = a * (lat-phi0s2+phi0s4-phi0s6); + + om0 = (double)1.0 - e2; + x = pow(om0,(double)0.5); + E1 = ((double)1.0 - x) / ((double)1.0 + x); + E2 = E1*E1; + E3 = E1*E2; + E4 = E1*E3; + A0 = (double)3.*E1/(double)2.-(double)27.*E3/(double)32.; + A1 = (double)21.*E2/(double)16.-(double)55.*E4/(double)32.; + A2 = (double)151.*E3/(double)96.; + A3 = (double)1097.*E4/(double)512.; + + + + tol = (double)1.e-5; + + dx = E - E0; + dy = N - N0; + M1 = AM0 + dy; + mu = M1 / (a*c0); + mus2 = A0 * sin((double)2.*mu); + mus4 = A1 * sin((double)4.*mu); + mus6 = A2 * sin((double)6.*mu); + mus8 = A3 * sin((double)8.*mu); + phi1 = mu + mus2 + mus4 + mus6 + mus8; + + if((((po2-tol)po2) + *phi=po2; + else if(*phi<-po2) + *phi=-po2; + + if(*lambda>GPS_PI) + *lambda -= p2; + if(*lambda<-GPS_PI) + *lambda += p2; + + if(*lambda>GPS_PI) + *lambda=GPS_PI; + else if(*lambda<-GPS_PI) + *lambda=-GPS_PI; + } + + *lambda = GPS_Math_Rad_To_Deg(*lambda); + *phi = GPS_Math_Rad_To_Deg(*phi); + + return; +} + + + + /* @func int32 GPS_Math_WGS84_To_ICS_EN ****************************** ** ** Convert WGS84 latitude and longitude to @@ -750,24 +1073,23 @@ void GPS_Math_Swiss_EN_To_WGS84(double E, double N, double *lat, double *lon) int32 GPS_Math_WGS84_To_ICS_EN(double lat, double lon, double *E, double *N) { - const double phi0 = 31.734090; - const double lambda0 = 35.212060; - const double E0 = 170251.0; - const double N0 = 1126868.0; - double phi, lambda, alt, a, b; + double const phi0 = 31.73409694444; // 31 44 2.749 + double const lambda0 = 35.21208055556; // 35 12 43.49 + double const E0 = 170251.555; + double const N0 = 1126867.909; + double phi, lambda, alt, a, b; - if (lat < 29.333) return 0; - if (lat > 33.28) return 0; - if (lon < 34.18) return 0; - if (lon > 37.67) return 0; - - a = GPS_Ellipse[27].a; - b = a - (a / GPS_Ellipse[27].invf); + int32 datum = GPS_Lookup_Datum_Index("Palestine 1923"); + int32 ellipse = GPS_Datum[datum].ellipse; + + a = GPS_Ellipse[ellipse].a; + b = a - (a / GPS_Ellipse[ellipse].invf); - GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, 124); - GPS_Math_Swiss_LatLon_To_EN(phi, lambda, E, N, phi0, lambda0, E0, N0, a, b); + GPS_Math_WGS84_To_Known_Datum_M(lat, lon, 0, &phi, &lambda, &alt, datum); + GPS_Math_Cassini_LatLon_To_EN(phi, lambda, E, N, + phi0, lambda0, E0, N0, a, b); - return 1; + return 1; } @@ -785,23 +1107,24 @@ int32 GPS_Math_WGS84_To_ICS_EN(double lat, double lon, double *E, ************************************************************************/ void GPS_Math_ICS_EN_To_WGS84(double E, double N, double *lat, double *lon) { - const double phi0 = 31.734090; - const double lambda0 = 35.212060; - const double E0 = 170251.0; - const double N0 = 1126868.0; - double phi, lambda, alt, a, b; - - - a = GPS_Ellipse[27].a; - b = a - (a / GPS_Ellipse[27].invf); - - GPS_Math_Swiss_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0, E0, N0, a, b); - GPS_Math_Known_Datum_To_WGS84_M(phi, lambda, 0, lat, lon, &alt, 124); + double const phi0 = 31.73409694444; // 31 44 2.749 + double const lambda0 = 35.21208055556; // 35 12 43.49 + double const E0 = 170251.555; + double const N0 = 1126867.909; + double phi, lambda, alt, a, b; + int32 datum = GPS_Lookup_Datum_Index("Palestine 1923"); + int32 ellipse = GPS_Datum[datum].ellipse; + + a = GPS_Ellipse[ellipse].a; + b = a - (a / GPS_Ellipse[ellipse].invf); + + GPS_Math_Cassini_EN_To_LatLon(E, N, &phi, &lambda, phi0, lambda0, + E0, N0, a, b); + GPS_Math_Known_Datum_To_WGS84_M(phi, lambda, 0, lat, lon, &alt, datum); } - /* @func GPS_Math_EN_To_LatLon ************************************** ** ** Convert Eastings and Northings to latitude and longitude diff --git a/reference/route/naviguide.gpx b/reference/route/naviguide.gpx index 13d1afbc2..ef82813bf 100644 --- a/reference/route/naviguide.gpx +++ b/reference/route/naviguide.gpx @@ -6,87 +6,87 @@ xmlns="http://www.topografix.com/GPX/1/0" xsi:schemaLocation="http://www.topografix.com/GPX/1/0 http://www.topografix.com/GPX/1/0/gpx.xsd"> - + - + 0.000000 A001 תחילת מסלול תחילת מסלול - + 0.000000 A002 - + 0.000000 A003 - + 0.000000 A004 - + 0.000000 A005 - + 0.000000 A006 - + 0.000000 A007 כניסה לשטח כניסה לשטח - + 0.000000 A008 - + 0.000000 A009 - + 0.000000 A010 - + 0.000000 A011 - + 0.000000 A012 - + 0.000000 A013 - + 0.000000 A014 - + 0.000000 A015 צומת עם שביל ירוק צומת עם שביל ירוק - + 0.000000 A016 - + 0.000000 A017 - + 0.000000 A018 - + 0.000000 A019 מערה -- 2.30.2